capture log close sublog
log using "$projdir/log/7_tables.txt", name(sublog) text replace

*-------------------------------------------------------------------------------
* Description: Make tables for the main text
* Authors: Justin Wiltshire, Denis Sosinskiy, Carl McPherson
*
* Updated: January 29, 2024
*-------------------------------------------------------------------------------
clear

* Ensure necessary folder exists
qui capture mkdir "$projdir/out/tables"

*-------------------------------------------------------------------------------
* Table 1: Main results through e21
*-------------------------------------------------------------------------------
tempfile core
tempfile core2

* Loop over each sample of treated units in the table
forval b = 0/2 {
	local alltrunits_`b'
	local n = 722513
	local end_tm = 21
	local y "avg_wkly_wage"
	foreach z in California NY {
				
		* Get list of treated units
		local filelist: dir "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'" files "*.dta"
			
		* Add the variables if needed
		foreach f of local filelist {
			local g
			if strpos("`f'", "_ate") == 0 {
				local g = subinstr("`f'", "tr_", "", .)
				local g = subinstr("`g'", ".dta", "", .)
				if (`b' == 0 & !inlist(`g', 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 1 & !inlist(`g', 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 2 & !inlist(`g', 6001, 6013, 6041, 6055, 6075, 6081, 6085, 6095, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) {
					local alltrunits_`b' "`alltrunits_`b'', `g'"
				}
			}
		}
	}
	
	local alltrunits_`b' = subinstr("`alltrunits_`b''", ", ", "", 1)

	* Get the 2010 population weights for each treated unit in sample b
	use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

	* Generate unique county id and keep only one observation per unit with pop10
	qui gen cty_fips = 1000*statefips + countyfips
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10

	* Merge in the MW growth
	qui merge 1:1 cty_fips using "$projdir/dta/build/cln/cty_mw_growth_e21.dta", nogen norep keep(3)
	
	* Keep only the treated units in sample b
	qui keep if inlist(cty_fips, `alltrunits_`b'')

	* Get the 2010-population-weighted % change in MW for b
	qui collapse (mean) pct_mw_growth [aw=pop10]
	qui sum pct_mw_growth
	local mw_grow_`b' = r(mean)
	di r(mean)

	* Loop over outcomes of interest
	foreach v in avg_wkly_wage employment {
		
		* Load the data
		if `b' == 0 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/agg_att/bias_corrected/bias_corr_att_distn.dta", clear
		}
		if `b' == 1 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/nolocal_att/bias_corrected/bias_corr_att_distn.dta", clear
		}
		if `b' == 2 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/nobaynyc_att/bias_corrected/bias_corr_att_distn.dta", clear
		}
		
		* Keep 2019Q4 for the treated aggregate only
		qui keep if _placeboID == 0 & _tm == `end_tm'

		* Capture and put in a matrix
		mat m_`b'_`v' = [.],[.]
		if `b' != 0 {
			mat m_`b'_`v' = m_`b'_`v'\m_`b'_`v'
		}
		qui sum gap_bc
		local gap_`v' = r(mean)
		local elast = `gap_`v''/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`gap_`v''],[.]
		mat m_`b'_`v' = m_`b'_`v'\[`elast'],[.]
		qui sum LB_95
		local lb_elast = r(mean)/`mw_grow_`b''
		qui sum UB_95
		local ub_elast = r(mean)/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`lb_elast'],[`ub_elast']
		qui sum p
		local pval = r(mean)
		mat m_`b'_`v' = m_`b'_`v'\[`pval'],[.]
		qui sum _Se
		local se_`v' = r(mean)
	}
		
	* Combine into two columns (one for each outcome)
	mat m_`b' = m_`b'_avg_wkly_wage,m_`b'_employment
	
	* Calculate OWE
	local owe = `gap_employment'/`gap_avg_wkly_wage'
	
	* Calculate OWE standard error using the delta method, then calculate OWE 95% CI
	mat jacobian = [1/`gap_avg_wkly_wage']\[-`gap_employment'/(`gap_avg_wkly_wage'^2)]
	mat varmat = [(`se_employment'^2)],[0]\[0],[(`se_avg_wkly_wage'^2)]
	mat owe_var = jacobian'*varmat*jacobian
	mata: st_matrix("owe_se", sqrt(st_matrix("owe_var")))
	mat owe_lb = `owe' - 1.96*owe_se
	mat owe_ub = `owe' + 1.96*owe_se

	* Create another column for own-wage elasticities and associated CIs
	if `b' == 0 {
		mat m_`b'_3 = [.],[.]\[m_`b'[2,3]/m_`b'[2,1]],[.]\[.],[.]\owe_lb,owe_ub\[.],[.]
	}
	if `b' != 0 {
		mat m_`b'_3 = [.],[.]\[.],[.]\[m_`b'[3,3]/m_`b'[3,1]],[.]\[.],[.]\owe_lb,owe_ub\[.],[.]
	}
	mat m_`b' = m_`b',m_`b'_3
}

* Combine into one matrix
mat results = m_0\m_1\m_2

* Make results tables for each industry
frmttable using "$projdir/out/tables/table_1", replace tex fragment statmat(results) sdec(2) ///
ctitles("","Average Weekly","Employment","Own-wage"\"","Earnings","","Elasticity") ///
rtitles("\textbf{All Treated Counties}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\"\textbf{Excluding Counties with}"\"\textbf{Local Minimum Wages}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\"\textbf{Excluding Counties in}"\"\textbf{SF Bay Area and NYC}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value")

*-------------------------------------------------------------------------------
* Table 2 (top two panels): Fast food estimates by estimator (SC)
*-------------------------------------------------------------------------------
tempfile core
tempfile core2

* Loop over each sample of treated units in the table
forval b = 0/1 {
	local alltrunits_`b'
	local n = 722513
	local end_tm = 21
	local y "employment"
	foreach z in California NY {
				
		* Get list of treated units
		local filelist: dir "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'" files "*.dta"
			
		* Add the variables if needed
		foreach f of local filelist {
			local g
			if strpos("`f'", "_ate") == 0 {
				local g = subinstr("`f'", "tr_", "", .)
				local g = subinstr("`g'", ".dta", "", .)
				if (`b' == 0 & !inlist(`g', 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 1 & !inlist(`g', 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) {
					local alltrunits_`b' "`alltrunits_`b'', `g'"
				}
			}
		}
	}
	
	local alltrunits_`b' = subinstr("`alltrunits_`b''", ", ", "", 1)

	* Get the 2010 population weights for each treated unit in sample b
	use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

	* Generate unique county id and keep only one observation per unit with pop10
	qui gen cty_fips = 1000*statefips + countyfips
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10

	* Merge in the MW growth
	qui merge 1:1 cty_fips using "$projdir/dta/build/cln/cty_mw_growth_e21.dta", nogen norep keep(3)
	
	* Keep only the treated units in sample b
	qui keep if inlist(cty_fips, `alltrunits_`b'')

	* Get the 2010-population-weighted % change in MW for b
	qui collapse (mean) pct_mw_growth [aw=pop10]
	qui sum pct_mw_growth
	local mw_grow_`b' = r(mean)
	di r(mean)

	* Loop over estimators of interest
	forval estimator = 0/1 {
		local estm bias_corrected
		local estm2 bias_corr_
		local estm3 bias_corrected
		local gap_name gap_bc
		if `estimator' == 1 {
			local estm uncorr
			local estm2 ""
			local estm3 ""
			local gap_name gap
		} 
		
		* Load the data
		if `b' == 0 {
			qui use "$projdir/dta/analysis/qcew/employment/naics`n'/agg_att/`estm3'/`estm2'att_distn.dta", clear
		}
		if `b' == 1 {
			qui use "$projdir/dta/analysis/qcew/employment/naics`n'/nolocal_att/`estm3'/`estm2'att_distn.dta", clear
		}
		
		* Keep 2019Q4 for the treated aggregate only
		qui keep if _placeboID == 0 & _tm == `end_tm'

		* Capture and put in a matrix
		mat m_`estm'_`b' = [.],[.]
		*if `b' != 0 {
			*mat m_`estm'_`b' = m_`estm'_`b'\m_`estm'_`b'
		*}

		qui sum `gap_name'
		local gap_`estm' = r(mean)
		local elast = `gap_`estm''/`mw_grow_`b''
		mat m_`estm'_`b' = m_`estm'_`b'\[`gap_`estm''],[.]
		mat m_`estm'_`b' = m_`estm'_`b'\[`elast'],[.]
		qui sum LB_95
		local lb_elast = r(mean)/`mw_grow_`b''
		qui sum UB_95
		local ub_elast = r(mean)/`mw_grow_`b''
		mat m_`estm'_`b' = m_`estm'_`b'\[`lb_elast'],[`ub_elast']
		qui sum p
		local pval = r(mean)
		mat m_`estm'_`b' = m_`estm'_`b'\[`pval'],[.]
		qui sum _Se
		local se_`estm' = r(mean)
	}
		
	* Combine into two rows(one for each estimator)
	mat m_`b' = m_bias_corrected_`b'\m_uncorr_`b'
}

* Combine into one matrix
mat results = m_0,m_1

* Make results table
frmttable using "$projdir/out/tables/table_2_top", replace tex fragment statmat(results) sdec(2) ///
ctitles("","Full","No Local MW"\"","Sample","Sample") ///
rtitles("\textit{A. Stacked synthetic control (Bias-corrected)}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\"\textit{B. Stacked synthetic control (uncorrected)}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value")

*-------------------------------------------------------------------------------
* Table 2 (bottom panel): Fast food estimates by estimator (DiD)
*-------------------------------------------------------------------------------
***************************************************
* Set sample 
***************************************************
clear all

qui capture mkdir "$projdir/dta/analysis/did"

gen outcome = ""
save "$projdir/dta/analysis/did/outcomes_did_samps.dta", replace 

use "$projdir/dta/build/cln/treated_cty.dta", clear 
gen treated = 1 

append using "$projdir/dta/build/cln/donor_cty.dta"

rename cty_fips area 

tempfile classify 
save `classify'

* Main data 
use "$projdir/dta/build/cln/analysis_panel_cty", clear

recode ind (722211=722513) if year <= 2011
recode ind (722110=722511) if year <= 2011
drop if employment == . 

sort area quarterly_date ind 
order area quarterly_date ind 

merge m:1 area using `classify', keep(3) nogen 
drop if statefip == 11 

***************************************************
* Prep for DiD 
***************************************************
clonevar cty_fips = area 
clonevar tm = quarterly_date 

keep if inrange(quarterly_date,199,251)

gen cs_tr = 0
replace cs_tr = 218 if statef == 6 //CA  
replace cs_tr = 216 if statef == 36 //NY

foreach v in "employment" "avg_wkly_wage" {
	gen temp = `v' if quarterly_date == 215
	egen temp2 = min(temp), by(area ind)
	gen rel_`v' = `v' / temp2 * 100 
	drop temp temp2 
}

forvalues l = 0/32 {	
	gen event_F`l' = (quarterly_date == cs_tr + `l')*treated 
}

forvalues l = 1/20 {
	gen event_L`l' = (quarterly_date ==  cs_tr -`l')*treated 
}

drop event_L1
order event???? 
drop quarter 

***************************************************
* CSDID 
***************************************************
foreach size in "all_large" "all_large_nolocal"  {
		local ind = 722513 
		foreach outcome in  "rel_employment" "rel_avg_wkly_wage" {
			
	preserve 
		
	di ""
	di " --- `size ' - `ind' - `outcome' --- "

	keep if ind == `ind'

	if("`outcome'" == "rel_avg_wkly_wage") local simple = "Average weekly earnings"
	if("`outcome'" == "rel_employment") local simple = "Employment"

	if("`outcome'" == "rel_avg_wkly_wage") local yrange = "-15(5)40"
	if("`outcome'" == "rel_employment") local yrange = "-30(5)20"

	* Geographies 
	if("`size'" == "all_large_nolocal") drop if inlist(area, 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005,36047,36061,36081,36085)

	* Divisors 
	if("`size'" == "all_large") local divisor = 58.53
	if("`size'" == "all_large_nolocal") local divisor = 49.31 

	keep if ind == `ind'

	di "" 
	di _N 
	di ""

	forv n = 1/5 {
	noi di ""	
	}
	noi di " CSDiD --- `size ' - `ind' - `outcome' --- "

	#delimit ;
	noi csdid `outcome' [w=pop10],
	 time(tm) 
	 ivar(cty_fips)
	 gvar(cs_tr) 
	 cluster(statefips) wboot rseed(7) method(dripw) agg(event)
	; 
	#delimit cr 

	mat rtable = r(table)'
	mat list rtable 

	xsvmat rtable, fast
	keep if inrange(_n,3,52)
	gen _tm = _n -17

	rename rtable1 beta 
	rename rtable5 lb 
	rename rtable6 ub 
	rename rtable2 stderr

	if("`size'" == "all_large") {
		di "`size' `outcome'"
		list beta lb ub if _n == _N 

		if("`outcome'" == "rel_avg_wkly_wage") {
			di "`size' `outcome'"
			local gap_avg_wkly_wage = beta[_N]
			local se_avg_wkly_wage = stderr[_N]
		}
		if("`outcome'" == "rel_employment") {
			di "`size' `outcome'"
			local gap_employment =  beta[_N]
			local se_employment = stderr[_N]		
		}
	}
 
	gen elast = beta/`divisor'
	gen ub_elast = ub/`divisor'
	gen lb_elast = lb/`divisor'
	gen samp = "`size'" 
	gen outcome = "`outcome'"

	keep _tm samp outcome beta *elast*
	order _tm outcome samp beta elast ub_elast lb_elast
	keep if inlist(_tm,21,33) 

	ds outcome samp _tm, not 
	foreach v in `r(varlist)' {
		replace `v' = round(`v',0.01)
	}

	gen ind = `ind'
	keep if _tm == 21
	drop _tm 
	gen time = 2019
	order time ind outcome samp beta elast ub_elast lb_elast
	sort time ind outcome samp 
	append using $projdir/dta/analysis/did/outcomes_did_samps
	save $projdir/dta/analysis/did/outcomes_did_samps, replace 

	restore 
	}

	* Calculate OWE
	local owe = `gap_employment'/`gap_avg_wkly_wage'

	* Calculate OWE standard error using the delta method, then calculate OWE 95% CI
	mat jacobian = [1/`gap_avg_wkly_wage']\[-`gap_employment'/(`gap_avg_wkly_wage'^2)]
	mat varmat = [(`se_employment'^2)],[0]\[0],[(`se_avg_wkly_wage'^2)]
	mat owe_var = jacobian'*varmat*jacobian
	mata: st_matrix("owe_se", sqrt(st_matrix("owe_var")))
	mat owe_lb = `owe' - 1.96*owe_se
	mat owe_ub = `owe' + 1.96*owe_se

	di "---`ind'---"
	di "`size' `outcome'"
	di "----------"
	di "`gap_employment' and `gap_avg_wkly_wage'"
	di "OWE: " round(`gap_employment'/`gap_avg_wkly_wage',0.01)
	di "LB: " 
	mat list owe_lb 
	di "UB: " 
	mat list owe_ub
	di "----------"
}

* Put estimates in a table
foreach s in "all_large" "all_large_nolocal" {		
		* Load the data
		qui use $projdir/dta/analysis/did/outcomes_did_samps.dta, clear
		
		* Keep sample of interest only
		qui keep if samp == "`s'" & outcome == "rel_employment"

		* Capture and put in a matrix
		mat m_`s' = [.],[.]

		qui sum beta
		local gap = r(mean)
		qui sum elast
		local elast = r(mean)
		mat m_`s' = m_`s'\[`gap'],[.]
		mat m_`s' = m_`s'\[`elast'],[.]
		qui sum lb_elast
		local lb_elast = r(mean)
		qui sum ub_elast
		local ub_elast = r(mean)
		mat m_`s' = m_`s'\[`lb_elast'],[`ub_elast']
	}

* Combine into one matrix
mat results = m_all_large,m_all_large_nolocal

* Make results table
frmttable using "$projdir/out/tables/table_2_bottom", replace tex fragment statmat(results) sdec(2) ///
ctitles("","Full","No Local MW"\"","Sample","Sample") ///
rtitles("\textit{C. Difference-in-differences}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} WBS CIs")

*-------------------------------------------------------------------------------
* Table 3 (top two panels): MCD main results through 2019
*-------------------------------------------------------------------------------
tempfile core
tempfile core2

* Loop over each sample of treated units in the table
forval b = 0/1 {
	local alltrunits_`b'
	local end_tm = 5
	local y "hwage"
	foreach z in California NY {
				
		* Get list of treated units
		local filelist: dir "$projdir/dta/analysis/mcd/`y'/`z'_subsample" files "*.dta"
			
		* Add the variables if needed
		foreach f of local filelist {
			local g
			if strpos("`f'", "_ate") == 0 {
				local g = subinstr("`f'", "tr_", "", .)
				local g = subinstr("`g'", ".dta", "", .)
				if (`b' == 0 & !inlist(`g', 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 1 & !inlist(`g', 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) {
					local alltrunits_`b' "`alltrunits_`b'', `g'"
				}
			}
		}
	}
	
	local alltrunits_`b' = subinstr("`alltrunits_`b''", ", ", "", 1)

	* Get the 2010 population weights for each treated unit in sample b
	use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

	* Generate unique county id and keep only one observation per unit with pop10
	qui gen cty_fips = 1000*statefips + countyfips
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10

	* Merge in the MW growth
	qui merge 1:1 cty_fips using "$projdir/dta/build/cln/cty_mw_growth_mcd_e21.dta", nogen norep keep(3)
	
	* Keep only the treated units in sample b
	qui keep if inlist(cty_fips, `alltrunits_`b'')

	* Get the 2010-population-weighted % change in MW for b
	qui collapse (mean) pct_mw_growth [aw=pop10]
	qui sum pct_mw_growth
	local mw_grow_`b' = r(mean)
	di r(mean)

	* Loop over outcomes of interest
	foreach v in hwage price {
		
		* Load the data
		if `b' == 0 {
			qui use "$projdir/dta/analysis/mcd/`v'/agg_att_subsample/agg_att_distn.dta", clear
		}
		if `b' == 1 {
			qui use "$projdir/dta/analysis/mcd/`v'/nolocal_att_subsample/stackedsc_att_distn.dta", clear
		}
		
		* Keep 2019Q4 for the treated aggregate only
		qui keep if _placeboID == 0 & _tm == `end_tm'

		* Capture and put in a matrix
		mat m_`b'_`v' = [.],[.]
		if `b' != 0 {
			mat m_`b'_`v' = m_`b'_`v'\m_`b'_`v'
		}
		qui sum p_gap
		local gap_`v' = r(mean)
		local elast = `gap_`v''/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`gap_`v''],[.]
		mat m_`b'_`v' = m_`b'_`v'\[`elast'],[.]
		qui sum LB_95
		local lb_elast = r(mean)/`mw_grow_`b''
		qui sum UB_95
		local ub_elast = r(mean)/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`lb_elast'],[`ub_elast']
		qui sum _Se
		local se_`v' = r(mean)
	}
		
	* Combine into two columns (one for each outcome)
	mat m_`b' = m_`b'_hwage,m_`b'_price
	
	* Calculate pass-through
	local pass = `gap_price'/(0.3*`gap_hwage')
	
	* Calculate pass-throug standard error using the delta method, then calculate OWE 95% CI
	mat jacobian = [1/(0.3*`gap_hwage')]\[-`gap_price'/(0.09*`gap_hwage'^2)]
	mat varmat = [(`se_price'^2)],[0]\[0],[(0.09*`se_hwage'^2)]
	mat pass_var = jacobian'*varmat*jacobian
	mata: st_matrix("pass_se", sqrt(st_matrix("pass_var")))
	mat pass_lb = `pass' - 1.96*pass_se
	mat pass_ub = `pass' + 1.96*pass_se

	* Create another column for pass-through and associated CIs
	if `b' == 0 {
		mat m_`b'_4 = [.],[.]\[`pass'],[.]\[.],[.]\pass_lb,pass_ub
	}
	if `b' == 1 {
		mat m_`b'_4 = [.],[.]\[.],[.]\[`pass'],[.]\[.],[.]\pass_lb,pass_ub
	}
	mat m_`b' = m_`b',m_`b'_4
}

* Combine into one matrix
mat results = m_0\m_1

* Make results table
frmttable using "$projdir/out/tables/table_3_top", replace tex fragment statmat(results) sdec(2) ///
ctitles("","Average","Price", "Pass-Through"\"","Hourly Wage","","") ///
rtitles("\textbf{All Treated Counties}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\textbf{Excluding Counties with}"\"\textbf{Local Minimum Wages}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs")

*-------------------------------------------------------------------------------
* Table 3 (bottom two panels): QWI Seperation Rates Results through e21
*-------------------------------------------------------------------------------
tempfile core
tempfile core2

* Loop over each sample of treated units in the table
forval b = 0/1 {
	local alltrunits_`b'
	local end_tm = 21
	local y sepr
	local n 7225
	local samp allage
	foreach z in California NY {
				
		* Get list of treated units
		local filelist: dir "$projdir/dta/analysis/qwi/`y'/naics`n'/`samp'/`z'" files "*.dta"
			
		* Add the variables if needed
		foreach f of local filelist {
			local g
			if strpos("`f'", "_ate") == 0 {
				local g = subinstr("`f'", "tr_", "", .)
				local g = subinstr("`g'", ".dta", "", .)
				if (`b' == 0 & !inlist(`g', 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 1 & !inlist(`g', 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) {
					local alltrunits_`b' "`alltrunits_`b'', `g'"
				}
			}
		}
	}
	
	local alltrunits_`b' = subinstr("`alltrunits_`b''", ", ", "", 1)

	* Get the 2010 population weights for each treated unit in sample b
	use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

	* Generate unique county id and keep only one observation per unit with pop10
	qui gen cty_fips = 1000*statefips + countyfips
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10

	* Merge in the MW growth
	qui merge 1:1 cty_fips using "$projdir/dta/build/cln/cty_mw_growth_e21.dta", nogen norep keep(3)
	
	* Keep only the treated units in sample b
	qui keep if inlist(cty_fips, `alltrunits_`b'')

	* Get the 2010-population-weighted % change in MW for b
	qui collapse (mean) pct_mw_growth [aw=pop10]
	qui sum pct_mw_growth
	local mw_grow_`b' = r(mean)
	di r(mean)

	* Loop over outcomes of interest
	foreach v in sepr sepnewr {
		
		* Load the data
		if `b' == 0 {
			qui use "$projdir/dta/analysis/qwi/`v'/naics`n'/`samp'/agg_att/bias_corrected/bias_corr_att_distn.dta", clear
		}
		if `b' == 1 {
			qui use "$projdir/dta/analysis/qwi/`v'/naics`n'/`samp'/nolocal_att/bias_corrected/bias_corr_att_distn.dta", clear
		}
		
		* Keep 2019Q4 for the treated aggregate only
		qui keep if _placeboID == 0 & _tm == `end_tm'

		* Capture and put in a matrix
		mat m_`b'_`v' = [.],[.]
		if `b' != 0 {
			mat m_`b'_`v' = m_`b'_`v'\m_`b'_`v'
		}
		qui sum gap_bc
		local gap_`v' = r(mean)
		local elast = `gap_`v''/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`gap_`v''],[.]
		mat m_`b'_`v' = m_`b'_`v'\[`elast'],[.]
		qui sum LB_95
		local lb_elast = r(mean)/`mw_grow_`b''
		qui sum UB_95
		local ub_elast = r(mean)/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`lb_elast'],[`ub_elast']
		qui sum p
		local pval = r(mean)
		mat m_`b'_`v' = m_`b'_`v'\[`pval'],[.]
		qui sum _Se
		local se_`v' = r(mean)
	}
		
	* Combine into two columns (one for each outcome)
	mat m_`b' = m_`b'_sepr,m_`b'_sepnewr
}

* Combine into one matrix
mat results = m_0\m_1

* Make results table
frmttable using "$projdir/out/tables/table_3_bottom", replace tex fragment statmat(results) sdec(2) ///
ctitles("","All Workers","Low-Tenure Workers") ///
rtitles("\textbf{All Treated Counties}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\ "\textbf{Excluding Counties with}"\"\textbf{Local Minimum Wages}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value")

*-------------------------------------------------------------------------------
* Table 4: Main results results through e33
*-------------------------------------------------------------------------------
tempfile core
tempfile core2

* Loop over each sample of treated units in the table
forval b = 0/2 {
	local alltrunits_`b'
	local n = 722513
	local end_tm = 33
	foreach y in avg_wkly_wage employment{
	foreach z in California NY {
				
		* Get list of treated units
		local filelist: dir "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'" files "*.dta"
			
		* Add the variables if needed
		foreach f of local filelist {
			local g
			if strpos("`f'", "_ate") == 0 {
				local g = subinstr("`f'", "tr_", "", .)
				local g = subinstr("`g'", ".dta", "", .)
				if (`b' == 0 & !inlist(`g', 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 1 & !inlist(`g', 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) ///
					| (`b' == 2 & !inlist(`g', 6001, 6013, 6041, 6055, 6075, 6081, 6085, 6095, 6097, 36005, 36047, 36061, 36081, 36085, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111)) {
					local alltrunits_`b' "`alltrunits_`b'', `g'"
				}
			}
		}
	}
	
	local alltrunits_`b' = subinstr("`alltrunits_`b''", ", ", "", 1)

	* Get the 2010 population weights for each treated unit in sample b
	use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

	* Generate unique county id and keep only one observation per unit with pop10
	qui gen cty_fips = 1000*statefips + countyfips
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10

	* Merge in the MW growth
	qui merge 1:1 cty_fips using "$projdir/dta/build/cln/cty_mw_growth_e33.dta", nogen norep keep(3)
	
	* Keep only the treated units in sample b
	qui keep if inlist(cty_fips, `alltrunits_`b'')

	* Get the 2010-population-weighted % change in MW for b
	qui collapse (mean) pct_mw_growth [aw=pop10]
	qui sum pct_mw_growth
	local mw_grow_`b' = r(mean)
	di r(mean)

	* Loop over outcomes of interest
	foreach v in avg_wkly_wage employment {
		
		* Load the data
		if `b' == 0 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/agg_att/covid_corrected/covid_corr_att_distn.dta", clear
		}
		if `b' == 1 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/nolocal_att/covid_corrected/covid_corr_att_distn.dta", clear
		}
		if `b' == 2 {
			qui use "$projdir/dta/analysis/qcew/`v'/naics`n'/nobaynyc_att/covid_corrected/covid_corr_att_distn.dta", clear
		}
		* Keep 2022Q4 for the treated aggregate only
		qui keep if _placeboID == 0 & _tm == `end_tm'

		* Capture and put in a matrix
		mat m_`b'_`v' = [.],[.]
		if `b' != 0 {
			mat m_`b'_`v' = m_`b'_`v'\m_`b'_`v'
		}
		qui sum p_gap
		local gap_`v' = r(mean)
		local elast = `gap_`v''/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`gap_`v''],[.]
		mat m_`b'_`v' = m_`b'_`v'\[`elast'],[.]
		qui sum LB_95
		local lb_elast = r(mean)/`mw_grow_`b''
		qui sum UB_95
		local ub_elast = r(mean)/`mw_grow_`b''
		mat m_`b'_`v' = m_`b'_`v'\[`lb_elast'],[`ub_elast']
		qui sum p
		local pval = r(mean)
		mat m_`b'_`v' = m_`b'_`v'\[`pval'],[.]
		qui sum _Se
		local se_`v' = r(mean)
	}
		
	* Combine into two columns (one for each outcome)
	mat m_`b' = m_`b'_avg_wkly_wage,m_`b'_employment
	
	* Calculate OWE
	local owe = `gap_employment'/`gap_avg_wkly_wage'
	
	* Calculate OWE standard error using the delta method, then calculate OWE 95% CI
	mat jacobian = [1/`gap_avg_wkly_wage']\[-`gap_employment'/(`gap_avg_wkly_wage'^2)]
	mat varmat = [(`se_employment'^2)],[0]\[0],[(`se_avg_wkly_wage'^2)]
	mat owe_var = jacobian'*varmat*jacobian
	mata: st_matrix("owe_se", sqrt(st_matrix("owe_var")))
	mat owe_lb = `owe' - 1.96*owe_se
	mat owe_ub = `owe' + 1.96*owe_se

	* Create another column for own-wage elasticities and associated CIs
	if `b' == 0 {
		mat m_`b'_3 = [.],[.]\[m_`b'[2,3]/m_`b'[2,1]],[.]\[.],[.]\owe_lb,owe_ub\[.],[.]
	}
	if `b' != 0 {
		mat m_`b'_3 = [.],[.]\[.],[.]\[m_`b'[3,3]/m_`b'[3,1]],[.]\[.],[.]\owe_lb,owe_ub\[.],[.]
	}
	mat m_`b' = m_`b',m_`b'_3
}
}

* Combine into one matrix
mat results = m_0\m_1\m_2

* Make results table
frmttable using "$projdir/out/tables/table_4", replace tex fragment statmat(results) sdec(2) ///
ctitles("","Average Weekly","Employment","Own-wage"\"","Earnings","","Elasticity") ///
rtitles("\textbf{All Treated Counties}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\"\textbf{Excluding Counties with}"\"\textbf{Local Minimum Wages}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value"\"\textbf{Excluding Counties in}"\"\textbf{SF Bay Area and NYC}"\"\hspace{2mm} Treatment Effect"\"\hspace{2mm} Elasticity"\"\hspace{2mm} Placebo-variance-based 95\% CIs"\ "\hspace{2mm} RMSPE-based \textit{p}-value")


* Close log file
log close sublog
